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Abstract. The development of consistent and stable quasicontinuum models for multi-dimensional 
crystalline solids remains a challenge. For example, proving stability of the force-based quasicontin- 
uum (QCF) model [s] remains an open problem. In ID and 2D, we show that by blending atomistic 
and Cauchy-Born continuum forces (instead of a sharp transition as in the QCF method) one ob- 
tains positive-definite blended force-based quasicontinuum (B-QCF) models. We establish sharp 
conditions on the required blending width. 



1. Introduction 



Atomistic-to-continuum coupling methods (a/c methods) have been proposed to increase the 
computational efficiency of atomistic computations involving the interaction between local crystal 
defects with long-range elastic fields (6}[7||l6|[l9|[2l|[25|_ 
such as the quasicontinuum model (denoted QCE f39] 



26,38 



Energy-based methods in this class, 
exhibit spurious interfacial forces ("ghost 



forces" ) even under uniform strain p|37| . The effect of the ghost force on the error in computing the 
deformation and the lattice stability by the QCE approximation has been analyzed in 
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|36||38]. 

An alternative approach to a/c coupling is the force-based quasicontinuum (QCF) approxima- 
tion (7l[l2l[l3[[23[[25], but the non-conservative and indefinite equilibrium equations make the iter- 



13 -15 . Indeed, it is an 



ative solution and the determination of lattice stability more challenging 
open problem whether the (sharp-interface) QCF method is stable in dimension greater than one. 
Many blended a/c coupling methods have been proposed in the literature, e.g., (I}|4| [T7|[22|[34| 



35 41 . In the present work, we formulate a blended force-based quasicontinuum (B-QCF) method, 
similar to the method proposed in f23^, which smoothly blends the forces of the atomistic and 
continuum model instead of the sharp transition in the QCF method. In ID and 2D, we establish 
sharp conditions under which a linearized B-QCF operator is positive definite. 

Firstly, we estab- 



23 



Our results have three advantages over the stability result proven in 
lish //^-stability (instead of iJ^-stability) which opens up the possibility to include defects in the 
analysis, along the lines of [15y30| . Secondly, our conditions for the positive definiteness of the lin- 
earized B-QCF operator are needed to ensure the convergence of several popular iterative solution 
methods for the B-QCF equations 14,24 . We note that the convergence of these popular iterative 



solution methods for the QCF equations cannot be guaranteed because of its indefinite linearized 
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operator [14||24| . Thirdly, our results admit much narrower blending regions, which is crucial for 
the computational efficiency of the method. 

The remainder of the paper is split into two sections: In Section [2] we analyze positivity of the 
B-QCF operator in a ID model, whereas in Section [3] we analyze a 2D model. Our methods and 
results are likely more widely applicable to other force-based model couplings. 

2. Analysis of the B-QCF Operator in ID 

2.1. Notation. We denote the scaled reference lattice by eZ := {ei : i G Z}. We apply a macro- 
scopic strain F > to the lattice, which yields 

yp :=FeZ = {F€e)eez. 
The space U of 2A^-periodic zero mean displacements u = {ug)g(^i from y_p is given by 



|u : ug+2N = Ui for ^ G Z, and I]^_7v+i ""^ = o|. 



and we thus admit deformations y from the space 

yp := {y : y = Yf + u for some u G U}. 

We set € = 1/N throughout so that the reference length of the computational cell remains fixed. 
We define the discrete differentiation operator, Du, on periodic displacements by 

[L)u)£ := , —CO < t < oo. 

We note that (Du)^ is also 2A^-periodic in i and satisfies the zero mean condition. We will denote 
{Du)^ by Due. We then define (D^^^u)^ and (D^^^u)^ for — oo < £ < oo by 

(2) (2) 

^(2) \ Due+i-Due _ (jj{3)\ _^ Dn^^Du^ 
/ e' e ' \ J e ' e 

To make the formulas more concise we sometimes denote Due by u'^, D^'^'^ue by n", etc., when there 
is no confusion in the expressions. 

For a displacement u € W and its discrete derivatives, we employ the weighted discrete £^ and 
i"^ norms by 



N \ 1/2 

|u||,2:=(e Yl l^^l'l ' ^=_^?\f,<^l^^l' 



\ e=-N+i / 
and the weighted inner product 

N 

(u,w) := ^ euewe- 
e=-N+i 

We will frequently use the following summation by parts identity: 

Lemma 2.1 (Summation by parts). Suppose {fk}^=m ^'^^ {9k}^=m ^'"^ sequences, then 

n m 

fk {gk+l - 9k) = [fn+ign+1 - fmOm] - ^ Qk+l {fk+1 - fk) ■ 
k=m k=n 

Also for future reference, we state a discrete Poincare inequality [3l] , 

||v||^oo < ||Z)v||£i for allv gU. 
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2.2. The next-nearest neighbor atomistic model and local QC approximation. We con- 
sider a one-dimensional (ID) atomistic chain with periodicity 2N , denoted y € 3^. The total 
atomistic energy per period of y is given by <f"(y) — ^Y1^=-n+i f^Vti where 

N 

i=-N+l 

for a scaled Lennard- Jones type potential (I8p8 and external forces fi. The equilibrium equations 
are given by the force balance at each atom: + = where 

We assume that the displacement u'^ = y° — yF is "small" and hence linearize the atomistic 
equilibrium equations about y^ to obtain 

(L'^u"), = /,, for £=-iV + l,...,iV, 

where {L°"v) for a displacement v G W is given by 

{L v)^ := (pp -2 + 02F ^2 • 

Here and throughout we use the notation (pp := (j)"{F) and (j)2p := 0"(2F), where (p is the potential 
in (2.1 ). We assume that (pp > 0, which holds for typical pair potentials such as the Lennard- Jones 
potential under physically relevant deformations. 

We will later require the following characterisation of the stability of L". 

Lemma 2.2. L*^ is positive definite, uniformly for N G N, if and only if cq := min(0^, </>^+4(/>2^) > 
0. Moreover, 

(L"u,u) > copull^a VuG^. 



Proof. The case (j)2p < was treated in 11 , hence suppose that 4)2p > 0. The coercivity estimate 
is trivial in this case, and it remains to show that it is also sharp. To that end, we note that 

Hence, testing with u'^ = (—1)^ (this is admissible since there is an even number of atoms per 
period), the second-neighbor terms drop out and we obtain (L"u, u) = 0^||Du|||2- D 



The local QC approximation (QCL) uses the Cauchy-Born extrapolation rule 38 ,39], that is, 
approximating + y'i_i in \2.1\ by 2y^ in our context. Thus, the QCL energy is given by 

N 

£''\y) = e + (2.3) 

e=-N+i 

We can similarly obtain the linearized QCL equilibrium equations about the uniform deformation 

'iqd^qd\ ^ e = -N + 1,...,N, 



where the expression of (L'^^'v)^ with v G W is 
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2.3. The Blended QCF Operator. The blended QCF (B-QCF) operator is obtained through 
smooth blending of the atomistic and local QC models. Let /3 : M — t- M be a "smooth" and 2-periodic 
blending function, then we define 

F^'^^iy) := miy) + (1 - MFr'iy), 

where F^'^^ is defined analogously to and := P{Fe£). Linearisation about yp yields the 
linearized B-QCF operator 

(L^^^^v), := f3e{L'^v)e + (1 - /3,)(L^^'v),. 

In order to obtain a practical atomistic-to-continuum coupling scheme, we would also need to 
coarsen the continuum region by choosing a coarser finite element mesh. In the present work we 
focus exclusively on the stability of the B-QCF operator, which is a necessary ingredient in any 
subsequent analysis of the B-QCF method. 

2.4. Positive-Definiteness of the B-QCF Operator. We begin by writing L*"?^-^ in the form 
L^i^f = (P'J.lI"^^ + ^'ipL\^^^ where 

j^i^'^^r] =e~'^ {-ve+i + 2vi - Vi-i) , and 



(2.4) 



(l^^'^^v)^ =/3^e-2 {-ve+2 + 2ve - ve-2) + (1 - M'^e-^ {-Vi+i + '^ve - ve-i) . 

Lemma 2.3. For any u & U, the nearest neighbor and next-nearest neighbor interaction operator 
can be written in the form 

(lJ'^'^-^u, u) =||Z?u||^2, and 

(L^-^^Vu) =[4||Z)u||| - e2||v/^Z)(2)^|||] + g + t, 
where the terms R and S are given by 

N N 

R= ^ 2e3Z)(2)/3, (z)n,)^ S= ^ e'^D^'^ P^D^^K^Due 
e=-N+i e=-N+i 

N 

and T= {D^^'^(3e+i)ueDue+i. 
e=-N+i 



(2.5) 



Proof. Since the proof of the first identity of Lemma 2.3 is not difficult, we only prove the identity 
for The main tool used here is the summation by parts formula. We note that 

(L2^-'u,u)= 2^ e/3i 2 ue + e{l-f3e) ^ ' 



-N+l 
N 



4 {-ue+i + 2ue - ui-i) 



■ — — -2 

l=-N+l 

N 



sr^ (-U£+2 + - Qui + 4u£_i - ui_2) 

+ 2^ ^2 

N 

=4pu|||+ e^pi(^-D^^^ue+i + D^^^Ui^ui. (2.6) 

£=-N+l 
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We then apply the summation by parts formula to the second term of (2.6) to obtain 

N 

e=-N+i 

N N 

e=-N+i e=~N+i 

We use the summation by parts formula again and change the index according to the periodicity 
so that we get 

N 



£=-N+l 



e=-N+i e=-N+i 

N N 



t=~N+l l=-N+l 
N 

e=-N+i 

N 

e=-N+i 

N 

e=-N+i 



(2.7) 



We now focus on the second term of (2.7). We repeatedly use the summation by parts formula to 
obtain 

AT 

l=-N+l 

N 

= Yl -^^Dp,[{Duif -{Dui^iY 

l=-N+l 

N 

N N 

e=-N+i e=-N+i 

N N N 

e=-N+i e=-N+i e=-N+i 

= R + S + T, 
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where R, S and T are defined in (2.5). 



Combining all of the above equalities, we obtain (2.4). 



□ 



We shall see below that the first group in (2.4) does not negatively affect the stability of the 
B-QCF operator. By contrast, the three terms R, S, T should be considered "error terms". We 
estimate them in the next lemma. 

In order to proceed with the analysis we define 

I:= {i£Z:0< Pi+j < 1 for some j G {±1, ±2}}, 
so that D(j)/3^ = for all i G {-N + 1, . . . iV} \ X and j G {1, 2, 3}, and K := jtX. 



Lemma 2.4. Let R, S and T be defined by (2.5), then we have the following estimates: 

IRI < e2||L>(2)/3| 



Du 



^2) 



and 



|S| < 2e2||D(2)/3||^^||z)u|||, 
|T| < e2^(Ke)^/2||^(3)^||^^ ii^^i 



(2. 



Proof. The estimate for R follows directly from Holder's inequality. 
To estimate S recall that D^'^^u^ := ^ which implies 



Therefore, S is bounded by 



N 



=-N+l 



|i?(2)u 



|L'u||^2 < 2e2||D(2)/3||£oo||Du| 



Finally, we estimate T by 
ITI = 



N 

E - 

i=-N+l 



< e 



'\\D^''^P\U^\\u\U2^j)\\Du\U2, 



We then apply the Holder inequality, the Poincare inequality and Jensen's inequality successively 
to ||u||£2(2:) to get 



'2 < i^e||Z?u||^i < 2Ke\\Du\\%. 



Therefore, we have 

ITI < e^IlD^^) 



Combining the above estimates, we have proven the second inequality in (2.8). 



□ 



We see from the previous result that smoothness of /? crucially enters the estimates on the error 
terms R, S, T. Before we state our main result in ID we show how quasi-optimal blending functions 
can be constructed to minimize these terms, which will require us to introduce the blending width 
into the analysis. The estimate (2.9) is stated for a single connected interface region, however, an 
analogous result holds if the interface has connected components with comparable width. A similar 
result can also be found in [19 . 
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Lemma 2.5. (i) Suppose that the blending region is connected, that is I = {1, . . . , K} without 
loss of generality, then (3 can be chosen such that 

WD^^^^Weo. <Cfs{Ker', for j = 1,2,3, (2.9) 

where Cp is independent of K and e. 

(ii) This estimate is sharp in sense that, if attains both the values and 1, then 

||Z)(^')/3||,oo > {Key^, for j = 1, 2, 3. (2.10) 

(iii) Suppose that J = {!,..., n} C X such that (3[V) = 0, (3{n) = 1 (or vice-versa), and 
0, n + 1 ^X, and suppose moreover that (2.9) is satisfied, then 

#{^ G J : D^'^Pi < -\{€Kr^] > ^K. (2.11) 

Proof, (i) The upper bound follows by fixing a reference blending function B £ C'^(R), B = 
in (-00,0] and S = 1 in [l,+oo), and defining /3{x) = B{{x - 2e)/{eK')) for K' = K - A. Then 
X = { 1 , . . . , K} , and a scaling argument immediately gives ( |2.9| ) . 

(ii) To prove the lower bound, suppose < Pi < 1 for £ = 1, . . . , Kq — 1, and /3o = and = 1- 
Then e'^f^^ j3[ = 1, from which infer the existence of Ki G {1, . . . , Kq} such that > l/(eii'o). 
This establishes the lower bound for j = 1. To prove it for j = 2 we note that, since = 1; 
P'ko+i — 0' ^^'^ hence we obtain 

e E /5" = /3xo+i-/3k<-V(e^o). 

We deduce that there exists K2 such that < -l/{e'^Ko{Ko - Ki)) < -l/{eKf. This implies 
(2.10) for j = 2. We can argue similarly to obtain the result for j = 3. 

(iii) Finally, to establish (2.11), let m G N be chosen minimally such that < —{eK)~'^ and 
/3o = 0; then m < n and we have 

1 ^ a" a" - flw ^ ^^'^/^ ^("^ ~ ^) 



>C-P'^ = eJ2(^T> 



where k := G : /3"' < —^{eK) Rearranging the inequality, we obtain 

1 > 1 ^ e{m-k) ^ ekCfs ^ kC^ 



2{tKY - {eKf 2{eKf " {eKf " K{eKY' 
and we immediately deduce that k/K > l/{2Cp), which concludes the proof of item (iii). □ 

We can summarize the previous estimates and get the following optimal condition for the size K 
of the blending region provided that (3 is chosen in a quasi-optimal way. Formally, the result states 
that L^'i'^f is positive definite if and only if ^ e^^/^. In particular, we conclude that the B-QCF 
operator is positive definite for fairly moderate blending widths. 



Theorem 2.1. LetX and K be defined as in Lemma 2.5. and suppose that j3 is chosen to satisfy 
the upper bound (|2.9[). Then there exists a constant Ci = Ci(C/3), such that 



{L'"i'^u,u) > (co - Ci\(t>'^p\ [K-'>/h-'/^])\\Du\\l Vu G U, (2.12) 
where cq = mm{(j)p, cpp + 4(/)2^) is the atomistic stability constant of Lemma 2.2. 
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Moreover, if Pi takes both the values and 1, then there exist constants 02,0^ > 0, independent 
of I, N, (f)"p and 4>2p, such that 



inf (L^^^/u, u) < 0'; + C2\cl>U - Cs\cl>'U [K-^/h~'/'] 



(2.13) 



ll^3u||^ 



Remark 2.1. Estimates (2.12) and (2.13) establish the asymptotic optimality of the blending width 
K 7z (-^1^ in the hmit as e 0: ( |2.12[ ) imphes that, if cq > and K > e"^/^, then L*"?^^ is 
coercive, while (2.13) shows that, if X ^ g^i/s then is necessarily indefinite. □ 

Proof. We first prove the lower bound. The blended force-based operator satisfies L^'^'^^ 

{L'i'fu,u) = A^||Du||| - e2,/.'2V||v^£'(2)u||| + (P'^pC^ + S + T) 

where Ap := (pp + 



b'^p. From Lemma 



2.4 



we have 



IR + S + TI < 



4||Z)(2)/3||,^ + (i^e)V2p(3)^||^^^ 



Since ||-D'^'''^/3||^c« < CfiiKt) ^ , so we have 



4iKe)-^ + {KeY'''{Ke 



1/2, 



|R + S + T| < Ce^ 
where we used the fact that K'"^ < K^^/'^e~^/'^ . 



\-3 



\Dm\ 



^-5/2g-l/2 



\Dvi\ 



£2) 



If (boTp < 0, then we obtain 



\Du 



|2 

1^2 . 



If (t)2p > 0, then 



{lI^'^u, u) = Ap\\Du\\l - e^^'ip\\^/pD^^^u\\l + (l)'^p{K + S + T) 



> 



i" 
hp 



■^-5/2^-1/2' 



\D\i\ 



12, 



(2.14) 



which is the corresponding result. 

To prove the opposite bound, let J be defined as in Lemma 2.5 (iii). We can assume this 
without loss of generality upon possibly shifting and inverting the blending function. We define 
J' := {£ £ J : D^^'^pi < -^{eKy^} and L := e#J' = aeK for some a > 1/(20,3), and a test 
function ^r through vq = ^ and 

and extending v'^ outside of I in such a way that ||L'v||^2 is bounded uniformly in I and N, and 
such that V is 2A^-periodic (see [l3] for details of this construction). 
With these definitions we obtain 

N 

T = e3 ^ Z)(3)/3,^,Z)^;,^, = e3 ^ ^(3)^^_^^/„^_^ 
f=-N+i eej' 



< 



4(ei^)3 

Recall that, by contrast, we have 



4ei^3 4 



^-5/2^-1/2^ 



IR + SI < C2K-^\\D^f\ 
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(a) Neighbor Set (b) Domain Decomposition 

Figure 1. (a) The 12 neighboring bonds of each atom, (b) The atomistic region is 
^la = Hex(ei?a). The blending region is Qj, = Eex{eRb) \ i^a- Here, Ra = 3, Rh = 7 
and K = 4. 



Combining these estimates, and using the fact that ||L)v||£2 is bounded independently of I and N, 



yields (2.13). □ 



3. Positive-Definiteness OF THE B-QCF Operator in 2D 

3.1. The triangular lattice. For some integer N ^ N and e := 1/A^, we define the scaled 2D 
triangular lattice 

r 1 1/2 

L := AgZ^, where Aq := [01,02] := e ^ 

where Oj, i = 1,2 are the scaled lattice vectors. Throughout our analysis, we use the following 
definition of the periodic reference cell 

n:=AQ{-N,Nf and C:=hnn. 

We furthermore set 03 = (— l/2e, \/3/2e)^, 04 := —01,05 := — 02 and og := —03; then the set of 
nearest-neighbor directions is given by 

TVi := {±01, ±02, ±03}. 

The set of next nearest-neighbor directions is given by 

7V2 := {±61, ±62, i&s}, where 6i:=oi+02, 62 := 02 + ^3 and 63 = 03-01. 

We use the notation M := A/i U M2 to denote the directions of the neighboring bonds in the 
interaction range of each atom (see Figure [T]). 

We identify all lattice functions v : L — t- with their continuous, piece affine interpolants with 
respect to the canonical triangulation T of with nodes L. 

3.2. The atomistic, continuum and blending regions. Let Hex(r) denote the closed hexagon 
centered at the origin, with sides aligned with the lattice directions 01,02,03, and diameter 2r. 

For Ra < Rb ^ N, we define the atomistic, blending and continuum regions, respectively, as 

:= Hex(e-Ra)) ■= Hex(ei?fc) \ Qa, and Qc '■= clos {Q \ {fla U i^b)) ■ 
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We denote the blending width by K := Rf, — Ra- Moreover, we define the corresponding lattice 
sites 

£":=£nl^a, C^:=Cnnh, and £^:=£nJlc- 

For simplicity, we will again use C as the finite element nodes, that is, every atom is a repatom. 
For a map v : L — t- and bond directions r, s G J\f, we define the finite difference operators 

v{x + r)-v{x) r> n ( \ Dsv{x + r) - Dsv{x) 

DrV[x) := and DrL>sV[x) := . 

e e 

We define the space of all admissible displacements, as all discrete functions L — )• which 
are Q-periodic and satisfy the mean zero condition on the computational domain: 

:= |u : L — )• : u{x) is fi-periodic and J2xi^c'^(^) ~ o}- 

For a given matrix B £ M^^^, det(i?) > 0, we admit deformations y from the space 

yB := {y : L — : y{x) = Bx + 'u(x), Vx G Lfor some u G W}. 

For a displacement u £U and its discrete directional derivatives, we employ the weighted discrete 
i'^ and i'^ norms given by 

/ Xl/2 

||u||^2 := I |n(a;)p I , ||u||£oo := max and 



/ 3 \ 1/2 

x£C i=l 



The inner product associated with is 

(u, w) := u{x) ■ w{x). 

3.3. The B-QCF operator. The total scaled atomistic energy for a periodic computational cell 
0, is 

^'^(y) =f E E ^(Drvix)) = 6^ E E [^(D^M^)) + </)(Z),,y(x))] , (3.1) 

x&C r&Af x&C i=l 

where (j) G C^(M^), for the sake of simplicity. Typically, one assumes (j){r) = ip{\r\); the more general 
form we use gives rise to a simplified notation; see also fso]. We define (p'{r) G and (/>"(r) G M^^^ 
to be, respectively, the gradient and hessian of (p. 

The equilibrium equations are given by the force balance at each atom, 

F%x;y) + fix;y) = 0, for x G £, (3.2) 

where f{x]y) are the external forces and F"'{x;y) are the atomistic forces (per unit volume e^) 

1 d£^{y) 



F^{x;y) :-- 



dy{x) 

3 , 3 



i=l i=l 
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Again, since u = y — y^, where ysix) = Bx, is assumed to be small we can linearize the atomistic 
equilibrium equation (3.2) about y^: 

{L''u''){x) = f{x), for xeC, 

where (L"v) (x), for a displacement v, is given by 

3 3 

(LV) (x) = - ^ ^"{Bai)Da,DaMx " " (j)" {Bh)Di,^Di,^v{x - bi), for x e C. 

i=l 1=1 

The QCL approximation uses the Cauchy-Born extrapolation rule to approximate the nonlocal 
atomistic model by a local continuum model 25 , 37 , 39 . According to the bond density lemma [30| 
Lemma 3.2] (see also [36j), we can write the total QCL energy as a sum of the bond density integrals 

where dry{x) = ^y{x + tr)\t=o denotes the directional derivative. We compute the continuum force 
F^{x;y) = — dy(x) ' ^^"^ linearize the force equation about the uniform deformation y^ to obtain 

(L^u^) (x) = /(x), for X e £. 

To formulate the B-QCF method, let the blending function : — )• [0,1] be a "smooth", 
0-periodic function. We shall suppose throughout that Ra, Rb are chosen in such a way that 

suvv{Da,,Da^^Da,.^P)clVLb Vi G {1, . . . , 6}^. (3.4) 

Then, the (nonlinear) B-QCF forces are given by 

F^'?^^(x; y) := /3(x)F'^(x; y) + (1 - l3{x))F%x; y), 

and linearizing the equilibrium equation F^'Jc/ _|_ j _ q about ys yields 

{L^^''fu^'i''f){x) = f{x), forxe£, 

where (L^'?^^v)(x) = /3(x)(LV)(x) + (1 - /3(x))(L'^v)(x). ^^'^^ 

Since the nearest neighbor terms in the atomistic and the QCL models are the same, we will 
focus on the second- neighbor interactions. We rewrite the operator L^'^'^^ in the form 

(L''«'=^v)(x) = ^(L^«'=^v)(x), 

where L^'?^^v(x) = /3(x)(L»(x) + (1 - /3(x))(L^v)(x), 
where the nearest-neighbor operators are given by 

La,^{x) = Ll.v{x) = -(j)"{Baj)Da^Da^v{x - Oj), 
and the second- neighbor operators, stated for convenience only for 6i = ai + 02, by 
(L^^u) (x) = - (j)"iBbi)Db,DbM^ - 61), while 
(Lg^u) (x) = - (t)"{Bhi) [Da^Da^u{x - ai) + Da^Da^u^x - 02) 

+ Da^Da2u{x - ai) + Da^Da2u{x - 02)] . 
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3.4. Auxiliary results. The following is the 2D counterpart of the summation by parts formula. 
The proof is straightforward. 

Lemma 3.1 (Summation by parts). For any u £h( and any direction r £ Z^, we have 

DrDru(x — r) ■ u{x) = — Dru{x — r) ■ Dru{x — r). (3.6) 

The second auxiliary result we require is a trace- or Poincare-type inequality to bound ||u||^2(Qj^) 
in terms of global norms. As a first step we establish a continuous version of the inequality we 
are seeking. The key technical ingredient in its proof is a sharp trace inequality, which is stated in 
Section O 

Lemma 3.2. Let ra < ri, ^ (0, 1/2], and let H := Hex(rb) \ Hex(ra); then there exists a constant C 
that is independent ofra,rt, such that 

\\u\\l2^jj-^<C[{rb-ra)rb\\ogrb\]\\du\\l2^^) £ H^Q), j^udx = 0. (3.7) 

Proof. Let S := 9Hex(l), and let dS denote the surface measure, then 

f-rb r 

I^lli2(/i') = ' ' W^^dSdr. 



Applying (5.1) with r^ = r and ri = 1 to each surface integral, we obtain 

WAhi^H) ^ i^b - ra){Co\\u\\l2^^) + Ci\\du\\l2^^)) , 



where Cq < Sr^ and Ci = 2rb\ logrbj. An application of Poincare's inequality yields (3.7). □ 



In our analysis, we require a result as (3.7) for discrete norms. We establish this next, using 
straightforward norm-equivalence arguments. 

Lemma 3.3. Suppose that < A/2, then 



u 



|2 ^ ry tr<o.,b^2 



< C {C''/Y\\DvL\\j2 VugU. (3. 



where C is a generic constant, and Cp := ^{eK)(eRb)\log{eR, 

Proof. Recall the identification of u with its corresponding Pi-interpolant. Let T G T with corners 
Xj, j = 1, 2, 3, then 

/ udx=^- — ^-'S^u(xi), and hence / udx = Vu G ZY. 
Jt ^ J^i 

Let ra ■= eRa and rj, := eRb, then H defined in Lemma 3.2 is identical to Qb- For any element 
T C fib it is straightforward to show that 

||u||^2(T) < C||n||i2(y). 

This immediately implies 

llull^fC^") < C\\u\\l2(^h), (3-9) 
for a constant C that is independent of e, Ra, K and u. Applying (3.7) yields 

||u|||(^f,) < C[{rb - ra)rb\ logrfe|] \\du\\l2^^y 
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Fix T £ T and let xj G T such that xj + aj £ T as well. Employing [30| Eq. (2.1)] we obtain 



3 3 

and summing over T £ T,T C (l we obtain that ||9ii||^2(f^') < C||Du||£2. This concludes the 
proof. □ 

3.5. Bounds on L^j^^^ • We focus only on the 6i-bonds, however, by symmetry analogous results 
hold for all second-neighbor bonds. As in the ID case, we begin by converting the quadratic 
form (L^^'^'^u, u) into divergence form. To that end it is convenient to define the bond-dependent 
symmetric bilinear forms and quadratic forms (although we write them like a norm they are typically 
indefinite) 

{r,s)b ■= r'^4>"{Bb)s, and \r\l := {r,r)h, for r, s, 6 E M^. 
Lemma 3.4. For any displacement u £U, we have 

{Llfu, u) = (Lg^u, u)-e^Yl - a2)\Da,DaMx - ai - 03)!^^ + Kb, + S;,,, (3.10) 

where 



Kbi := - {Dail3{x - 2ai){Da^u{x - 2ai), Da^DaiUix - oi - 02)) 



61 



+ Da^Pix - a2){Daiu{x - ai),Daj^Da2u{x - ai - a2))^J, and (3.11) 
Sfei := - Dai -0^1/3(2; - 2ai)(u(x - ai), Da^Da^uix - ai - 02))^,^. 



xeC 



Proof. For this purely algebraic proof we may assume without loss of generality that (f)"{Bbi) = I. 
In general, one may simply replace all Euclidean inner products with (•, 



Starting from (3.5), we have 

{L^fu, u) ={Llu, u) + (L^^u - L^^u, /3u) 

=m^u, u) - ^ P{x)u{x) ■ [Db^Db^u{x - bi) - Da^Da^u{x - ai) 

-Da^Da^uix - 02) - DaiDa2u{x - Ol) - Da^Da^uix - 02)] . 

We will focus our analysis on (i^^u — Lg^u, /3u). 

Noting that 61 = oi + 02, one can recast Db^Db^u{x — bi) as 

Db^Db^u{x - bi) 

= 4t [uix + 61) - 2u{x) + u(x - 61)] 

=DaiDa2u{x) + Da^Da^u{x - ai) + Da^Da^uix - 02) + Da^Da^uix - Ql - 02). 
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Applying the summation by parts formula (3.6) to (-^^^^u — L'^_^u, /3u), we get 

(L^iU - Ll^U, f3u) = - ^ /3(x)u(x) • Da^Da-^Da^U^X - ai) - Da^Da^Da2u{x - ai - 02 

x&C 

= - e"^ ^ /3(x)u(x) • Da^Da^Da^Da^uix - Oi - 02) 

=e^ ^ Da-^Da.,Da2u{x - oi - 02) • Dai (/^(^^ " ai)u{x - ai)^ 

=e^ ^ Da-^Da^Da^uix - ai - 02) • /3(x)i:>aiii(2; - ai) + u{x - ai)Da^l3{x - ai) 



xec 



Another application of the summation by parts formula ( |3.6[ ) converts (-^^^^u — Lg^u,/3u) into 
- Ll^u,(3u) =e'^'^Da^Da2Da^u{x - ai - 02) • {u{x - ai)Da^l3{x - ai)) 



x&C 



^ Da^Da2u{x - ai - 02) • {Da^Pix - a2)Daiu{x - ai)) 
- ^ Da^Da^uix - ai - 02) • - a2)Da^Da2u{x - oi - 02)) . 



The first two terms on the right-hand side can be rewritten as 

X] [Da^Da^Da^uix - ai - 02) • (n(x - ai)Da^l3{x - aij) 



x&C 



- Da^Da2u{x - ai - 02) • {Da^f^ix - a2)Da^u{x - ai))} 
-e^ ^ (^(a^ - ai)Da^Daj^(3{x - 2ai) j • Da^Da^uix - oi - 02) 



e"^ ^ {Z)ai/3(a:; - 2ai)Da^u{x - 2ai) ■ Da2Da2u{x - ai - 02) 



+ Da^Pix - a2)Da^u{x - ai) ■ Daj^Da-Mx - ai - 02)} 



Thus, we obtain (3.10) and (3.11) 



□ 



Next, we will bound the singular terms R;,^ and S^^, for which we introduce the notation 



max \\Da,Da^f3\Uoo, and p^^^/^H^. 

l<i,j<6 



max \\Da^Da^DaJ 
l<i,j,k<6 " » J 



Lemma 3.5. The terms H-bi o-nd S^j defined in (3.11) are bounded by 



and 



a, 6 



(3.12) 
(3.13) 



where C is a generic constant and Cp^ is defined in Lemma 



5.< 



Proof. According to the expression of Rf,^ given in (3.11) and noting that 



2 4 ..o 1,2 4 ..n 

\Da2Da2u\\^2 < -^\\Du\\^2 and \\Daj^Da2u\\^2 < -^\\Du\\^2, 



POSITIVE-DEFINITENESS OF THE BLENDED FORCE-BASED QUASICONTINUUM METHOD 



15 



we immediately obtain the first inequality of (3.12). 
We first rewrite Sf,^ as 



Sfti = - X] Da^Da^f3{x - 2ai){Da^Da^u{x - ai - 02), u{x - ai))^^ 
xec 

= - e'^'^Da^Da^f3{x - 2ai)Da2{Da2u{x - ai - a2),u{x - ai - 02))^,^ 

+ e'' ^ Da^Da^l3{x - 2ai){Da.,u{x - ai - 02), Da^uix - oi - 02))^^ 
=e'^'^ Da^Daj^Da^Pix - 2ai - a2){Da2u{x - ai -02) ■ u{x - ai -02))^^ 
e^^^Da^Da^P{x - 2ai)\Da2u{x -ai- 02) 



|2 



(3.14) 



x<^C 



For the second term in (3.14), we have 

Y Da,DaJ{x - 2ai)\DaMx -ai- a2)\l^ | < e^W {Bhi)\\\D^^^ pu| 

For the first term, we have 

^e^Y^Da2Da^Da^l3{x - 2ai - a2){Da2u{x - ai - a2),u{x - ai - 02))^ 

< e2 \<l)"{Bhi)\ \\MDa,Da,DaMil\\DAp^ < W'{Bhi)\ p(3)/3||,^ WAi^iC) 



The last inequality comes from the assumption (3.4), which ensures that supp(Da2^ai-C'a^/3) C 
Applying Lemma 3.3 yields the bound for Sb^. □ 

To summarize the estimates of this section we define a self-adjoint operator L by 

3 

(Lu, u) := {L^vl.u) - e^^^/5(a; - a2)\Da^Da^+^u{x - ai - 02)]^^; (3.15) 



then. Lemma 3.4 and Lemma |3.5| immediately yield the following result. 



Corollary 3.1. Suppose that Ra and are defined such that (3.4) holds; then, for all u ^U, 
(L^^^^u,u) > (Zu,u) -C7C"[e2||D/3||^oo +e2||D(2)/3||^^ +e2^^'^||^(3)^||^^jp^||2^^ (3_^g) 



where C is a generic constant, C" := maxj=i_2,3 \ (t)"{Bhj)\ and Cp^ is defined in Lemma 

Based on the analysis and numerical experiments in [30] for a similar linearized operator we 
expect that the region of stability for L is the same as for that is, L is positive definite for a 
macroscopic strain B if and only if L"" is positive definite. However, we are at this point unable to 
prove this result. Instead, we have the following weaker result. The proof is elementary. 

Proposition 3.1. Suppose that B G M^^^ is such that is positive definite, 



5.< 



(L^u,u) > 7c||Du||| VuGZ^, 
and suppose that (j)"{Bbj) < 61 where 6 < 7c/4, then L is positive definite, 



{Lu,u) >^\\Du\\% yueU, 



(3.17) 



with 7 = 7c — 4(5. 
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3.6. Positivity of the B-QCF operator in 2D. The blending width K is again a crucial ingredi- 
ent in the stability analysis for Lf't'^f . Due to the simple geometry we have chosen it straightforward 
to generalize Lemma 2.5 to the two-dimensional case, using the same arguments as in ID. 

Lemma 3.6. It is possible to choose P such that 

\\D^^'^Ph^<Cp{Ke)-K for j = 1,2,3, (3.18) 

Since we cannot fully characterize the stability of L in terms of the stability of or we 
will only prove stability of L^'?'^'^ subject to the assumption that L is stable. Proposition 3J gives 
sufficient conditions. 



Theorem 3.1. Suppose that (3 is chosen quasi- optimally so that (3.18) is attained; then, 

{L'"''fu,u)>-fbqcf\\Du\\J,, 

where 

li>,cf -.= 1-00" [e-'l'K-'l^\eR,\og{eR,)\'% 
where C is a generic constant and C" is defined in Corollary \3.1\ 

In particular, if L is positive definite (3.17) and if K is sufficiently large, then L^i'^^ is positive 
definite. 



Proof. From Corollary |3 . 1 1 and (3.18) we obtain 



(L^'^^Vu) > {^-CC" [e2(ei^)-i + e2(eK)-2 + e2(^^)-5/2|^^^l^g(^^^)|i/2]|||^^||| 

Remark 3.1. Suppose that 7 > 0, uniformly as — )• 00 (or, e — )• 0). In this limit, we would like 
to understand how to optimally scale K with Ra. (Note that Ra controls the modeling error; cf. 
Remark |3.3[ ) We consider three different scalings of Ra. 

Case 1: Suppose that Ra is bounded as e — )■ 0. In that case, we obtain 

Ibqcf - CC" e-'/^K-^^MRa + K) \og{e{Ra + K))\^'^ 

= -CC" (1 + f ) ( log(ei^) + log(l + f )) I'/' 

- - CC" K-^\\og{eK)\^l^. (3.19) 

From this it is easy to see that L'"?'^-^ will be positive definite provided we select K ^ | loge|^/^. 



Case 2: Suppose that 1 <^ Ra ^ e ^; to precise, let Ra 



~ e 



for some a G (0, 1). Then, a 



similar computation as (3.19) yields 

lb,cf - 7 ^ K-'^'^K + e-°)(log e + log(i^ + e"")) | 

and we deduce that, in this case, L^i^f will positive definite provided we select K S> e~"/^| loge|^/^. 

Case 3: Finally, the case when the atomistic region is macroscopic, i.e., Ra = 0{e~^), can be 
treated precisely as the ID case and hence we obtain that, if we select K ^> e~^/^, then l}"i^f is 
positive. 

In summary, we have shown that, in the limit as e — )• 0, if L is positive definite, Ra = 0(e~°) 
and if we choose 

loge|i/4. 



a 



0, 



loge|l/5g-a/5^ 
.-1/5 



< a < 1, 
a = 1, 



(3.20) 
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Figure 2. Visualization of the construction discussed in 3.2: the white region is 



the atomistic domain, the hght gray region the blending region, the medium gray 
region and dark gray regions together are the set J and the dark gray region is the 
set J'. 



0. We emphasize that 
□ 



then the B-QCF operator L^^c/ positive definite and ^^^^^ ~ 7 as e 
these are very mild restrictions on the blending width. 

It remains to show that the sufficient conditions we derived to guarantee positivity of are 
sharp. A result as general as (2.13) in ID would be very technical to obtain; instead, we offer a 
brief formal discussion for a special case. 

Remark 3.2. We consider again the limit as e — t- 0, and for simplicity restrict ourselves to the case 
where <C ~ and <C Ra ~ for O<0<Q!<1. In particular, R^ ~ e~" as well. 

We assume that Da^fiix) = for all x G J" C as depicted in Figure [2] The set J should be 
chosen so that its size is comparable with that of C'^ , but sufficiently small to still allow j3 to satisfy 
the bound (3.18). We can now repeat the ID argument along atomic layers to obtain that 

-3 _ 



-3+36 



Da,Da,Da,P{x) < -^{eKY 

for all 2; in a subset J' d J containing entire atomic planes, that has comparable size to J\ that 
is, # J' ~ KR^ ~ e"^-". 

Suppose now that (Bh\) has a negative eigenvalue A with corresponding normalized eigenvector 
u G M^, then we seek test functions of the form u{x) = fj,{x)u. It is now relatively straightforward, 
applying the ID argument in normal direction and using a smooth cut-off in the tangential direction, 
to construct fi supported in J' with Da2fJ-{x) ~ {e'^#J'')~^^'^ so that ||Du||^2 ~ 1, and 



^ Da2DajDa^P{x - 2ai - a2){Da2u{x - ai- 02), u(x - ai 



02) 



< 



e^XiY.Da^Da.DaJix 

xec 



2ai - a2)Da2n{x - ai - a2)n{x - ai - 02) 



-1/2 



,(50-q)/2 



This shows that, if X <C e~"/^, then L**?"^-^ is necessarily indefinite. 

In summary, for the specific interface geometry and a particular choice of (3 (which does, however, 
lead to the quasi-optimal bound (3.18)) we have shown that Theorem 3.1 is sharp up to logarithmic 
terms. □ 
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Remark 3.3. In practise, for the computation of different types of defects, we would first choose 
an appropriate scahng Ra = e~" for the atomistic region, considering the accuracy of the B-QCF 
method, and then choose the blending width K in order to ensure stability. 

For instance, for a point defect in 2D with zero Burger's vector it is expected that the displace- 
ment field satisfies Ua{x) = ya{x) — Bx ~ e/r, where r is the distance from the defect (30p3|. With- 
out coarse-graining, the local continuum (QCL) model has a modeling error of order 0{e^\d^Ua\) 



(see [lO 20,29] for proofs in ID and |40| for a proof in arbitrary dimensions); and although we 



have not established it rigorously, we expect that modeling error for the B-QCF method outside 



the atomistic region is also of second order; see also 13 . 

From u{x) e/r we can make the reasonable assumption that |3'^ya| ~ e/?^^, from which we 
obtain (assuming also stability) that the total error is of the order 

\\d{ya - ybqcf)\\L^ ~ e^\\d^ya\\mn\n,) ~ / r\r'^\'^drV''^ ~ R-^ . 

Hence, if we wish to obtain ||5(ya — yhqcj)\Lp- ~ e'^, < /c < 3, then we need to choose 

Ra ~ e"''/^, and consequently K > e~^^'^^\ loge|^/^. 

With this choice we can ensure both the stability and O(e^) accuracy of the B-QCF method; 
provided that our assumption that the B-QCF method has indeed a second-order modelling error 
is correct. □ 

4. Conclusion 

We have studied the stability a blended force-based quasicontinuum (B-QCF) method. In ID 
we were able to identify an asymptotically optimal condition on the width of the blending region 
to ensure that the linearized B-QCF operator is coercive if and only if the atomistic operator is 
coercive as well. In the 2D B-QCF model, we have obtained rigorous sufficient conditions and have 
presented a heuristic argument suggesting that they are sharp up to logarithmic terms. In 2D our 
proof of coercivity of relies on the coercivity of the auxiliary operator L defined in (3.15), for 
which we cannot give sharp conditions at this point. 

The main conclusion of this work is that the required blending width to ensure coercivity of the 
linearized B-QCF operator is surprisingly small. 

Our analysis in this paper is the first step towards a complete a priori error analysis of the 
B-QCF method, which will require a coercivity analysis of the B-QCF operator linearized about 
arbitrary states, as well as a consistency analysis in negative Sobolev norms. 

5. Appendix: A Trace Inequality 

In the following trace theorem, 5(1) denotes the unit sphere in M'', r := |x| and := x/|x|. Upon 
taking = 1 and employing standard orthogonal decompositions it is easy to check that the result 
is sharp. In particular, for d = 2, consider the case u{x) = log 

Lemma 5.1. Let d>2,ip: S{1) — )• (0,1] be Lipschitz continuous, and S := {ip{a)a : a G ^(l)}. 
Moreover, let < ro < ri < 1, and A := Uro<r<ri ('^^)' ^^^''^ 

ll^lli^CroE) < Co||n||i2(^) + Ci||an||i.(^), Vn e H\A), (5.1) 

where Cn = f— 1 , and Ci = < '^L ^ (5.2) 

n-roVn; \2ro/(d-2), d > 3. ^ ' 
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Proof. Since A is a Lipschitz domain we may assume, without loss of generality that u G C^{A). 
The symbol dS denotes the {d — l)-dimensional Hausdorff measure in W^. 
Let ro < s < ri, then 



\u\^dS = r^-^ / |n(rofT)|2d5^ 



u{sa) 



r=ro 



d^ 
dr 



u{ra)dr 



dSa 



<2r^-^ [ \u{sa)\'^dS^ + 2r^-^ [ [ du ■ adr 

Js Jt. Jr=ro 



dSa- 



(5.3) 



By hypothesis we have |c7| < 1 for all cr G S, hence the second term on the right-hand side can be 
further estimated, applying also the Cauchy-Schwartz inequality, by 

2 



2r: 



d-l 



du ■ adr 



r=ro 



dSa < 2r^~i 



E J r=ro 



r-'^+^dr / r'^-^\du{r(T)\'^drdS„ 



r=ro 



2r',-\j{s) - Jiro)) 



r=ro J rS 



\du\'^dSdr 



<2ro^-^(J(s)-J(ro))||an||i,(^), 

where J'{t) = t^^+i, that is, J{t) = logi if d = 2 and J{t) = t~'^+^/i-d + 2) if d > 3. Since J{s) 
is negative and strictly increasing for s < 1 we obtain 

2 



2r" 



du ■ adr 



r=ro 



(5.4) 



dS,<2r^,-'\Jiro)\\\du\\i,^^y 
Inserting (5.4) into (5.3), multiplying the resulting inequality by s*^^^ and integrating over s € 



(ro,ri) yields 



^f-^O |L,I|2 



L2(roS) 



s=ro 



< 2r\ 



„(i-i 



s=ro 



\u\ dS ds 



\u{sa)\^dS^ ds + 2r^-^J{ro) '^ 



Dividing through by ''^ we obtain 



li2(roE) 



Finally, estimating r^ — rf > (ri — ro)r^ ^ yields the stated trace inequality. 
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